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Abstract 

We propose an efficient numerical scheme for the resolution of a non-hydrostatic 
Saint-Venant type model. The model is a shallow water type approximation of the 
incompressbile Euler system with free surface and slightly differs from the Green- 
Naghdi model. 

The numerical approximation relies on a kinetic interpretation of the model and 
a projection-correction type scheme. The hyperbolic part of the system is approx¬ 
imated using a kinetic based hnite volume solver and the correction step implies 
to solve an elliptic problem involving the non-hydrostatic part of the pressure. 

We prove the numerical scheme satishes properties such as positivity, well¬ 
balancing and a fully discrete entropy inequality. The numerical scheme is con¬ 
fronted with various time-dependent analytical solutions. Notice that the numerical 
procedure remains stable when the water depth tends to zero. 

Keywords : shallow water flows, dispersive terms, finite volumes, finite differences, 

projection scheme, discrete entropy 
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1 Introduction 

Non-linear shallow water equations model the dynamics of a shallow, rotating layer of 
homogeneous incompressible fluid and are typically used to describe vertically averaged 
flows in two or three dimensional domains in terms of horizontal velocity and depth 
variations. The classical Saint-Venant system [8] with viscosity and friction [251 [26] l36] 
is particularly well-suited for the study and numerical simulations of a large class of 
geophysical phenomena such as rivers, lava flows, ice sheets, coastal domains, oceans 
or even run-off or avalanches when being modihed with adapted source terms [IHl [TTl 
[M] . But the Saint-Venant system is built on the hydrostatic assumption consisting in 
neglecting the vertical acceleration of the fluid. This assumption is valid for a large 
class of geophysical flows but is restrictive in various situations where the dispersive 
effects - such as those occuring in wave propagation - cannot be neglected. As an 
example, neglecting the vertical acceleration in granular flows or landslides leads to 
signihcantly overestimate the initial flow velocity [3^ [32] , with strong implication for 
hazard assessment. 

The derivation of shallow water type models including the non-hydrostatic effects 
has received an extensive coverage [28l [20l [9] [SS] [391 [13 [H] cind numerical techniques 
for the approximation of these models have been recently proposed [33 [23 [13 ED- 
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In [T7], some of the authors have presented an original derivation process of a 
non-hydrostatic shallow water-type model approximating the incompressible Euler and 
Navier-Stokes systems with free surface where the closure relations are obtained by a 
minimal energy constraint instead of an asymptotic expansion. The model slightly dif¬ 
fers from the well-known Green-Naghdi model j2^. The purpose of this paper is to 
propose a robust and efficient numerical scheme for the model described in |T^. The 
numerical procedure, based on a projection-correction strategy |22], is endowed with 
properties such as consistency, positivity, well-balancing and satishes a fully discrete 
entropy inequality. We emphasize that the scheme behaves well when the water depth 
tends to zero and hence is able to treat wet/dry interfaces. As far as the authors know, 
few numerical methods endowed with such stability properties have been proposed for 
such dispersive models extending the shallow water equations. 

The paper is organized as follows. In Section we recall the non-hydrostatic model 
proposed in im and we give a rewritting of the system. A kinetic description of the 
model is given in Section and it is used to derive the numerical procedure and to prove 
its properties that are detailed in Sections and In Section we first study the 
semi-discrete schemes (in space and in time) and then we establish some properties of 
the fully discrete scheme. In Section we prove the entropy inequality for the fully 
discrete scheme. Stationary/transient analytical solutions of the model are proposed in 
Sectionj^and hnally the numerical scheme is confronted with analytical and experimental 
measurements. 


2 A depth-averaged Euler system 

Several strategies are possible for the derivation of shallow water type models extending 
the Saint-Venant system. A usual process is to assume potential flows and an extensive 
literature exists concerning these models UDIEBESIEIEIEI]. An asymptotic expansion, 
going one step further than the classical Saint-Venant system is also possible IT8| H2] 
but such an approach does not always lead to properly dehned and/or unique closure 
relations. In this paper, we start from a non-hydrostatic model derived and studied 
in HU. where the closure relations are obtained by a minimal energy constraint. 

The non-hydrostatic model we intend to discretize in this paper has several interesting 
properties 

• the model formulation only involves first order partial derivatives and appears as 
a depth-averaged version of the Euler system, 

• the proposed model is similar to the well-known Green-Naghdi model but keeps a 
natural expression of the topography source term. 
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2.1 The model 


So we start from the system (see Fig. [^for the notations) 


OH d . . ^ 


dt 

d 


d 


r^2 . 9 


-(Hu) + — {Hu^ + = -{gH + 


dt dx 

d d 

^AHw) +—{Hum) = 2p^,^, 


dzh 


dx 


dt 


dx 


3(m) 


dx dx 

We consider this system for 


( 1 ) 

( 2 ) 

(3) 

(4) 


t > to and X G [0, L], 

u = {u, w)'^ denotes the velocity vector and the non-hydrostatic part of the pressure. 
The total pressure is given by 

P=|^ + Pnh- (5) 

The quantities (n, uJ, p) correspond to vertically averaged values of the variables (u, tc, p) 
arising in the incompressible Euler system. 



Figure 1: Notations: water depth H{x,t), free surface H + Zb{x,t) and bottom Zb{x,t). 


The smooth solutions H, u, w, p^h of the system 0-0 also satisfy the energy 
balance 

— (r; + gHzb) + ^ (u{p + gHzb + = 0, (6) 


where 


V = 


H{v? + w'^) g 


2 


(7) 
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We can rewrite 0 under the form 



dr] dG 

(8) 

with 

rj = r] + gHzb, 

G = u(rj + , G = G + gHzbU, 

G — G + (9) 

In the sequel, we 

will also use the dehnitions 



rjhyd = + 

(10) 


Vhyd Vhyd “1“ QH 

(11) 


Ghyd = u{r]hyd + ^H'^) , 

(12) 


Ghyd = u{fihyd + ^H'^) , 

(13) 


which are functions of the unknowns. 

The system Q-Q is completed with initial and boundary conditions that will be 
precised later. 


Remark 2.1 Notice that simple manipulations of Eqs. Q and (|^ lead to the relation 

'H^ + 2Hzb\ d rH^ + 2Hzb 


m 


A 

dx 


-u 


= Hw, 


(14) 


and Eq. Q could be replaced by {lA). In this paper we use (|^ which leads to keep 
the analogy with the divergence operator in the Navier-Stokes equations as shown in the 
following paragraph. 


2.2 A rewriting 

Let us formally define the operator Vgw by 


^SW f 


- 2 / 


(15) 


that is a shallow water version of the gradient operator. Likewise, we define a shallow 
water version of the divergence operator div^^ under the form 


d{Hu) d{H + 2zb),, 
diy„„ u = — - u --h 2w, 


dx 


dx 


(16) 


where u = (u, w)'^. In dehnitions ( [l^ and ( [l^ , we assume the considered quantities 
are smooth enough. The dehnition of the operators Vsw and div^^ implies we have the 
identity 


WswP-u dx = [Hup] 


ai 


pdivsi„ u dx, V p, Vu, 


(17) 
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where I is any interval of M. Notice that Vsw and divg^o are H and Zb dependent operators 
and when necessary we will use the notations Vsw (•; H) and div^^ (.; H). 

The system Q-(|^ can be rewritten under the compact form 


OH d ^ 


dt 

d 


^ iu.Hu) + Vo = 

div^^u = 0, 


-gHVoZb, 


(18) 

(19) 

( 20 ) 


with the notation 


Vo/ = 


Then the system (18 )-(20) appears as a Id shallow water version of the 2d incompressible 
Euler system. 


2.3 Pressure equation 


For H > 0, Eq. (19) can also be written in the nonconservative form 


5u _5u ^ rr 1 

^ ^ ^^ 0 ^^ + J^^swPnh = -g^oZb, 


( 21 ) 


and applying (16) to Eq. (21) leads together with (20) to the relation 


_ f^dp^h 
dx V dx 


1 L + , f8(H + 2z,)Y\ 

^nr-" dx^ + dx 




_2d'^Zb , + Zb) ^ dzbd{H + Zb) 


+ gH- 


dx^ 


2g 


dx dx 


( 22 ) 


Notice that Eq. (22) also reads 


/ \ ^ 

As^p^b = ‘2Hi^j +2u 


.-2^'^^b d‘^{HZb) dzbd{H + Zb) 


+ gH- 


dx'^ 


2g 


dx dx 


(23) 


with 


Agw divg^o I j^^sw 


Conversely, Eq. (21) with (23) give the divergence free condition (20). The resolution of 


Eq. (23) - requiring the inversion of a non local operator - gives the expression for the 


non-hydrostatic pressure term p^^- ^ discrete approximation of will be dehned in 


paragraph 4.3 and used for the numerical solution of (18)-(20). 
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Remark 2.2 In all the writings of the model, the pressure term p^y, appears as the 
Lagrange multiplier of the divergence free condition. As in the incompressible Euler 
system, it is not possible to derive a priori bounds for the pressure terms. And hence, it 
is possible to obtain nonpositive values for the total pressure p defined by (|^. 

Such a situation means the fluid is no longer in contact with the bottom and the 
formulation of the proposed model is no longer valid since the bottom of the fluid has to 
be considered as a free surface. 

Even if the proposed model can be modified to take into account these situations, we 
do not consider them in this paper and we will propose in paragraph \4 . 6| a modification 
allowing to ensure that the total pressure remains nonnegative. 


2.4 Other formulations 


One of the most popular models for the description of long, dispersive water waves is 
the Green-Naghdi model [2H]. Several derivations of the Green-Naghdi model have been 
proposed in the literature [2H1 EH SHI [37]. For the mathematical justihcation of the 
model, the reader can refer to ^m\ and for its numerical approximation to inminiEii 


Introducing a parameter a and starting from the system (18)-(20), we write 

A 

dx ' 


OH d ^ 

~A+' EL (■^“) “ 


dt 

^(Hu) + gj (u.Hu) + Vo ( ifC ) + VZp„k = -g-ffVozo. 
div,“„u = 0, 


(24) 

(25) 

(26) 


with 


and 


V“ f = 

^sw J 


Tjdf_ I d{H+2zp) 
^ dx^ dx 


f 


a d{Hu) 
div,,„ u = —- - u 


-af 
d{H + 2zb) 


+ aw. 


(27) 


dx dx 

The value a = 2 gives exactly the model Q-(|^. The system (24)-(26) is completed 
with the energy balance 


d^a 




where 


fj^ = E I 


2a-1 




+ lH^ + gHz,. 


Following [3lj (see also ini). the Green-Naghdi model with flat bottom reads 

dH d ^ 


(28) 

(29) 

(30) 
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d{Hu) 


(31) 

dt 

d 

) + = 2 ^ 3 -’ 

(32) 

d{Hu) 

dH 3 

(33) 

dx 


corresponding to (24)-(26) with the value a = 3/2. And hence it appears that the pro¬ 
posed model and the Green-Naghdi system only differ from the value of a in Eqs. (24 ),(26). 
Notice that the fundamental duality relation 


Pnhdivlu dx = [Hup^f^gi - / V," dx, 


holds for any interval I. 

It seems to the authors that the consistency of relations Eqs. (27) and (29) with 
the divergence free condition and the energy balance for the 2d Euler system is only 
obtained for the value a = 2 and in this paper we mainly focus on the choice a = 2 i.e. 
on the model (18)-(20). But the numerical procedure proposed and described in this 
paper is also valid for the model corresponding to any nonnegative value of a. 


3 Kinetic description 

In this section, we propose a kinetic interpretation for the system 0-0 completed 
with 0- The kinetic description will be used in Section to derive a stable, accurate 
and robust numerical scheme. 

The kinetic approach consists in using a description of the microscopic behavior 
of the system jlO]. In this method, a fictitious density of particles is introduced and 
the equations are considered at the microscopic scale, where no discontinuity occurs. 
The kinetic interpretation of a system allows its transformation into a family of linear 
transport equations, to which an upwinding discretization is naturally applicable. 

Following PI, we introduce a real function y defined on M, compactly supported 
and which has the following properties 

f x{-w) = x{w) > 0 

I = 1- ^ 

Among all the functions y satisfying ( [3^ , one plays an important role. Indeed, the 
choice 

1 / ^2x1/2 

= (35) 

with = max(0,a;), allows to ensure important stability properties [411 [5]. In the 
following, we keep this special choice for y. 
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3.1 Kinetic interpretation of the Saint-Venant system 

The classical Saint-Venant system m m\ corresponds to the hydrostatic part of the 
model , it reads 




completed with the entropy ineqnality 

dTjhyd ^ dGhyd 


< 0 , 


(36) 

(37) 

(38) 


dt dx 
with fihyd and Ghyd dehned by @.(1131. 

Let us construct the density of particles M{x,t,^) playing the role of a Maxwellian: 
the microscopic density of particles present at time t, at the abscissa x and with velocity 
^ given by 

1 / \ 1/2 
gn V 


M{H,u,0 = -X (— 

c \ c 


(39) 


with c = \ i G M. The equilibrium defined by (39) corresponds to the classical ki 


netic Maxwellian equilibrium, used in mi for example. It satishes the following moment 
relations, 


M{H,u,0d^ = 


H 

Hu 




eM{H,%Odi = Hu^ + g — . 


(40) 


The interest of the particular form ( [39| ) lies in its link with a kinetic entropy, see [5] 
where the properties of Hx{f,^,z) are studied, refers to the kinetic entropy used 
in [S]. Consider the kinetic entropy. 


HkU, z) = ^f + + gzf, 

where / > 0, ^ G M and z G M, and its version without topography 

7f/f.o(/.0 = k/ + Lrl"' 

2 D 

Then one can check the relations 


(41) 


(42) 


(43) 


C,HK(M(H,u,£,),C„Zt)dC, = G^ya. 


(44) 


These dehnitions allow us to obtain a kinetic representation of the Saint-Venant 
system HU- 


9 








Proposition 3.1 The pair of functions {H, Hu) is a strong solution of the Saint-Venant 

es the kir 

dzb dM 


system (36)-(37) if and only if M{H,u,^) satisfies the kinetic equation 

dM 


(B) 


e- 


dM 


9- 


= Q, 


dt dx dx d^ 
for some “collision term” Q{x,t,^) which satisfies, for a.e. {x,t), 


(45) 



iQdi = 0. 


(46) 


Proof of prop. 3.1| Using (40), the proof relies on a very obvious computation. 


Remark 3.2 The proposition 3.1 remains valid if instead of (35), the equilibrium M is 
built with any function satisfying (34). 


This proposition produces a very useful consequence. The non-linear shallow water 
system can be viewed as a single linear equation for a scalar function M depending 
nonlinearly on H and u, for which it is easier to find simple numerical schemes with 
good theoretical properties. 


3.2 Kinetic interpretation of the depth-averaged Euler system 

Since we take into account the non-hydrostatic effects of the pressure, the microscopic 
vertical velocity 7 of the particles has to be considered and we now construct the new 
density of particles M{x, t, 7 ) defined by a Gibbs equilibrium; the microscopic density 
of particles present at time t, abscissa x and with microscopic horizontal velocity ^ and 
microscopic vertical velocity 7 is given by 




(47) 


where 6 is the Dirac distribution and c = y 
Then we have the following proposition. 

Proposition 3.3 For a given p.^^, the functions {H,u,w) satisfying the divergence free 
condition (16), are strong solutions of the depth-averaged Euler system described in (j^- 
(§,(@ if and only if the equilibrium M{x,t,^,'y) is solution of the kinetic equations 




nh ) 


dM dM 

+ e- 


dt 


dx 


. ‘^Pnh\ dzb I d 


dM 

9M 
H d-f 


— Qnh, (48) 
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where Qnh = Qnh{x,t,^,'y) is a “collision term” satisfying 

I Qnh d^d'y / ^Qnh d^d'j / 'yQnhd^d'y 0 . 

J^2 J^2 

Additionally, the solution is an entropy solution if 

a + ^2 2 2 \ 

--1-+ Q^h ) Qnhd'ydf < 0. 


( 49 ) 


(50) 

Notice that in the case of the Saint-Venant system, the particular choice of M defined 


by (39) ensures 


, 2^2 


^ ^ 2 Qdi — 0 . 


Proof of prop. 3.3| From the definitions ([M)),( [47|) a nd the proof results from easy 
computations, namely by integrating the relation (48) 


/ (Bnh) d^d'y, / ^{Bnh) d^d'y, and / 'y{Bnh) d^d'y. 
7r2 7r2 7r2 

Likewise, the energy balance is obtained calculating the quantity 


e+1^ 


2 2 \ 

^-M^ + gzA {BnQ d^dj. 


Because of the special role of the equation (|^, it is not easy to describe it at the kinetic 
level. 


4 Numerical scheme 


In this section we propose a discretization for the system (18)-(20). In order to process 


step by step, we first establish some properties for the semi-discrete schemes in time and 
then in space. Then we study the fully discrete scheme. 


For the sake of simplicity, the notations with are dropped. We write the system (18 )- 


(19) in a condensed form 


with 


X = 


f)X f) 


S(.Y) = 




/ Hu 

Hu 

, F{X) = 

Hu^ + IH^ 

Hw / 


\ Huw 


0 

-gHVoZb 


(51) 


(52) 


and 


Fn.h 


0 

^sw Pnh 


divergence free condition (20) is satisfied. 


with VswPnh defined by (15). The expression of pnh is defined by (23) and ensures the 
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4.1 Fractional step scheme 

For the time discretization, we denote t”' = Ylk<n where the time steps At^ will 
be precised later though a CFL condition. Following j22], we use an operator splitting 
technique resulting in a two step scheme 

(53) 


At^ 

^n+l _ j^n+ 1/2 




KV = 0 - 


(54) 


The non-hydrostatic part of the pressure is defined by (23) and ensures, as already 
said, that the divergence free constraint (20) is satished i.e. 


diw,„ = 0 


(55) 


The discretization of Eq. (23) is given hereafter. The system (53)-(54) has to be com¬ 


pleted with suitable boundary conditions that will be precised later, see paragraph 4.3.3 


The prediction step (53) consists in the resolution of the Saint-Venant system and a 
transport equation for i.e. 


^n+l/2 


,d{Huy 


dx 


{Hu, 


,n+l/2 _ 


d 


= {Hu)^ -Ar—{ Hu^ + - At^gH^^, 


= {Hw)^ - At 

and the correction step (|54j) writes 


dx 
d{Hwu 


.dzh 


dx 


dx 


^rt+l _ j^n+1/2 


with 


U 


= 


Un+l/2 _ 

Ar 

Plf 

ffn+l 

(FFm)”+i 

{Hw)^+^ 

y 

ffn+l ’ 

Hn+l 

) • 


(56) 

(57) 

(58) 

(59) 

(60) 


More precisely, due to the expression of the operator Vsw given in (15), the notations 
means and the same remark holds for the operator div^^. 

Then inserting satisfying (60) in relation (55) gives the governing 


n+1 


equation for p"^ 


diy,, 


(. 




^ v.„p;;+d = d-diy. 


At/ 


(Hu 


jn+ 1/2 


(Hw 


|n+l /2 


^n+ 1/2 ’ ^n+ 1/2 


(61) 


that is a discrete version of (23). Notice also that in Eq. (60) we have used the fact that 

^n+l ^ ^n+l/2^ 

appears that the right hand side of (|61[) can be evaluated by the 
conservative variables {H, Hu, given by (56)-(58) , the hrst step of the time 

scheme. 
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Proposition 4.1 The scheme (53)-(55) satisfies a semi-discrete (in time) entropy in¬ 
equality of the form 


d 


~n+l ^ ^ | \\F'(X^) + S{X\ , 


with 


Tin , \ n 

fj- = fj{X^) = — + (w-)2 j + |(i7-)2 + gH-Zb, 


defined by (§, 


Proof of prop. |4.1| Multiplying Eq. (57) by , we obtain after classical manipulations 

-n+l/2 _~n _ f ^ |(^”)')) + 

fjn+1/2 

+ (62) 


Vhyd ^hyd 


with 


Hn 


Vlyi = + Iw"? + 

defined by (0. Likewise, multiplying Eq. ([5^ by leads to 


j^n+l/2 


[w 


jEfn f) / jEfn 


jjn+l/2 


[W 


.n+l/2 _ 


w^Y- (63) 


Notice that the last term appearing in Eq. (62) and in Eq. (63) is non negative. These 
error terms are due to the explicit time scheme. The sum of the two previous equations 
gives the inequality 


~n+l/2 <gn_ ^ (ArYOi ||F'(X") + S(X 

ox 


”)ll 2 


(64) 


Now we multiply (60) by (i7u)”+^ and after simple manipulations it comes 

Rn+l 


TTn+1/2 / a 

(m»+ 1)2 = - Ai" ( ^((i/«)"+VS‘) 


+fi‘ {Yjhu) 


d 


n+l _ ^n+l rn'^+l ^ 2Zb) 

OX 


jjn+1/2 


(65) 


and 

Rn+l 


[w 


n+l\2 


Tjn+l/2 0 - 71 + 1/2 

+ 2 Ai>::++"-^‘ - ( 


( 66 ) 
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Using the two previous equations and (55) gives the inequality 


In Eqs. (|65|) and (|66|), the error terms due to the time discretization are non-positive. 

(67) 


d 


.n+l < ~n+l/2 _ ^_ 


Finally Eq. (67) coupled with Eq. (64) gives the result. 


4.2 The semi-discrete (in space) scheme 


To approximate the solution X = {H, Hu, Hw)"’" of the system (51), we use a finite 
volume framework. We assume that the computational domain is discretized with I 
nodes Xi, i = 1,... ,1. We denote Cj the cell of length Axi = Xj+ 1/2 — Xj_i /2 with 
Xi+ 1/2 = {xi + Xi+i)/2. We denote Xi = {Hi, q,^^i, q^^iY with 




Ax. 


X{x, t)dx. 


i JCi 


the approximate solution at time t on the cell Ci with q^^i = HiUi, qz,i = HiWi. Likewise 
for the topography, we define 


^b,i 


1 

Axi 



Zh{x)dx. 


The non-hydrostatic part of the pressure is discretized on a staggered grid (in fact the 
dual mesh if we consider the 2d case) 


r^i + l 


Pnh,i-\-lf2 


Axi 


i+1/2 J Xi 


Pnh{x,t)dx, 


XXj_|_i/2 Xi- 

Now we propose and study the semi-discrete (in space) scheme approximating the 
model (51) and the divergence free condition ( [20| ). The semi-discrete scheme writes 

dX, 


Axi 


dt 


+ (Xi+l/2- — T)-1/2+) + Rnh,i — 0, 


dW^,i+l/2 ({Uj}) = 0, 


( 68 ) 

(69) 


where (69) is a discretized version of the divergence free condition (20) which we detail 
below and with the numerical fluxes 


75+1/2+ — X(w , Xj+ 1 , Zb^i, ^b,i+l) + iSj+i/ 2 + 

Xi+1/2- = X(W ) Xj+1, Zb^i, Zb^i-\-l') T iSj+i/2—• 

X is a numerical flux for the conservative part of the system, 5 is a convenient dis¬ 


cretization of the topography source term, see paragraph 4.3 
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Since the first two lines of (53) correspond to the classical Saint-Venant system, the 
nnmerical fluxes 

/ FH,i+l/2 \ 

Fi+i/2± = ( Fg^^i^i/2± I , (70) 

V Fq^^i+i/2 ) 

can be constructed using any numerical solver for the Saint-Venant system. More pre¬ 
cisely for FH-,i+i/2,77j^,i+i/2± we adopt numerical fluxes suitable for the Saint-Venant 
system with topography. Notice that from the dehnition (52), since only the second com¬ 
ponent of S{X) is non zero, only Fq^ has two interface values under the form i^q^,j+i/ 2 ±- 
For the dehnition of F)j^,i+i/ 2 , the formula (see jB]) 


Fq^,i+l/2 — FH^i+l/2Wi+i/2, 


with 


Wi+l/2 — 


(71) 

(72) 


Wi if FH,i+ii2 > 0 
Wi+i if FH,i+ii2 < 0 

can be used. 

Combining the hnite volume approach for the hyperbolic part with a hnite difference 
strategy for the parabolic part, the non-hydrostatic part Rnh,i is dehned by 


Rrih.i 


0 

^sw,i Pnh 


where the two components of Vsw,i Pnh are dehned by 
XXi Xsw,i Pnh 


XXi XsWjiPnhlo 


FI-i{pnh,i+l/2 Pnh,i—ll2) 

(Ci+l Ci) F 1/2 (Ci Ci—l) 5 

^XXi-^-i/2Pnh,i+l/2 F XXi—X/2Pnh,i—l/2^^ 


(73) 

(74) 


with 


0 = 


Hi “ 1 “ 


And in (69), diy 5 i„_j+i /2 (u) is dehned by 

Aa;i+i/ 2 div^^,i+i /2 (u) = {Hu)i+i - {Hu)i - (n* -1- Mi+i)(0+i - Q) 

+ Axi+i/2 {wi+i + Wi) . (75) 


Notice that in the dehnitions (73)-(74) and in the sequel, the quantity pnh means {pnh,j}- 
Likewise in Eq. (75) and in the sequel, u means {u^}. 

In a hrst step, we assume we have for the resolution of the hyperbolic part i.e. 
the calculus of F)+i/ 2 -,F)_i/ 2 +, a robust and efficient numerical scheme. Since this step 
mainly consists in the resolution of the Saint-Venant equations there exists several solvers 
endowed with such properties (HLL, Rusanov, relaxation, kinetic,...), see |12) . 

We assume we have for the prediction step a numerical scheme which is 
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(i) consistent with the Saint-Venant system (36)-(37), 
(n) well-balanced i.e. at rest = X'^ in (53), 

(in) satisfying an in-cell entropy of the form 


Xx. 


dTjhyd,^ 

dt 


+ [Ghyd,i+l/2 — Ghyd,1-1/2] < 0, 


with 


H, 


Vhyd,i = + f + gHiZb, 


and Ghyd,i+i /2 is the entropy flux associated with the chosen hnite volume solver. 
Then the following proposition holds. 


Proposition 4.2 The numerical scheme (68),(69) 

(i) is consistent with the model 0 -@> 

(a) preserves the same steady state as the lake at rest, 

(Hi) satisfies an in-cell entropy inequality associated with the entropy fi{t) analogous to 
the continuous one defined in 


(i^i+i /2 ~ ^ 1 - 1 / 2 ) < di, in Gi, 


(76) 


with 


Hi Vi.^i) Vhyd,i T Hi ^ ; 

Gi+l/2 = Ghyd,i+l/2 + -^n',i+l/ 2 'M^i+l/ 2/2 + {Hu)i+i/2Pnh,i+l/2- 
and di is an error term satisfying di = 0{Ax^), 

(iv) ensures a decrease of the total energy under the form 

d 


dt 


AzjT/i < 0. 


(77) 


The inequality (76) is obtained by multiplying (scalar product) the two momenta equa¬ 
tions of (68) by Uj which corresponds to a piecewise constant discretization. For the 
hyperbolic part it allows to derive a semi-discrete entropy, see nig. 

The error term di in the r.h.s. of (76) comes from the discretization of the non¬ 
hydrostatic part corresponding to the incompressible part of the model. In order to 
eliminate di, a more acurate discretization of the velocity held - in accordance with the 
approximation of the divergence free condition - would be necessary. 
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Proof of prop. 4.2| (i) Since we have assumed that the numerical scheme for the pre¬ 
diction part is consistent, we have 

7{X,X,z,z) = F{X). 

Likewise S is a consistent discretization of the topography source term. It is easy to 
prove the non-hydrostatic terms given by (73),(74) is a consistent discretization of Rnh 
that proves the result. 

(a) When Uj = (0,0)^ for j = i — l,i,i -\- 1, the hyperbolic part being discretized 
using a well-balanced scheme we have 

Fi+l/2- = -Pi-l/2+ = 0, 


dX, 

dt 


= 0 , 


and the scheme (68),(69) reduces to 

Rnh,i = (0,0,0)'^, 

ensuring the scheme is well-balanced. 

(Hi) Multiplying the first two eguations of system (68) by the first two components of 
fj'iXi) with 




gH. 


Ui 

Wi 


we obtain 


Axi 


^Vhydj'i 

dt 


W ( G^/i 3 /d,i+l/ 2 — 1/2+j 4“ ^sw,i Pnh\i A 0. 


(78) 


In Eg. (78), the three first terms are obtained as in W- The proof of theorem 2.1 in 
can be used without any change, except for the vertical kinetic energy, namely 

2 ’ 

that is not considered in since the model is hydrostatic. In order to obtain the con¬ 
tribution of the vertical kinetic energy in Eg. (78), we proceed as follows. 

Multiplying the third component of eguation of (68) by Wt i.e. the third component 
offjfiXi) we obtain 


Axi 


dHiWi 

dt 


4“ ( RH,i+l/2^i+l/2 l/2^i—1/2 4“ AXi Xsm^iPnfi\2^i 0- 


The first term in the preceding eguation also writes 


dHiWi 

dt 


d (H, 

Wi = ^ I —U!i I + 


wf dHi 

~2~dt 


(79) 
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For the fluxes, it comes 


{FH,i+l/2Wi+l/2 — — FH,i+l/2 


^1+1/2 


Fh i- 


H,i-l/2- 


+ FH,i+l/2Wi+i/2 (wi - _ FH,i-l/2Wi-i/2 (wi - • (§ 0 ) 


^h/2 

2 

Wi 


Using the first equation of (68) and the definition (72), the sum of Eqs. (79) and (80) 
gives 


d fH,J\ , ^ ^ ^ _ 

Q-j- ( 2 2+1/2 2 1/2 2 /^Xi \s'w,iPnh\2^i 

,*+1/2] _ ('U^i+l ~ Wi)'^ — - [FH^i-1/2 

with the notations [a]+ = max{a,0), [a]_ = min{a,0) a = [a]+ + [a]_. Therefore it 
comes 


]^{Wi- Wi-if, 



, d [Hi fi\ ^ w, 


i+1/2 


Fn.i- 


W 


H,i-ll2- 


i-1/2 




and the left hand side of the preceding equation is exactly the contribution of the vertical 
kinetic energy over the energy balance (78). 


Adding (78) to (81) gives 


+ ( G^j+ 1 / 2 - ~ Gi-i/2+ ) + 


Hi 

Wi 


•^sw,i Pnh ^ 0 , 


(82) 


with GiJ^i/ 2 - = Ghyd,i-\-ii 2 - + -Fk, 2 + 1/2 ci'rid it remains to rewrite the last terms in 

Eq. (82). Using the definitions (73);(74)j we have 


^Xi ^swfiPnh\2^i 
^Xi ^sw,i 


with 


^Ax2-i-i/2Pn/7, 2 + 1/2 “ 1 “ -^^ 2 —l/ 2 Pn/i, 2 — 1 / 2 ^ ^25 (^3) 

dii(^Pnh,i-^l /2 Pn/2,2—1/2)^2 

“l“Pr 2 / 2 , 2 +l /2 (C2+I 0)^2 ^22/2,2— 1/2 (0 0 —1)^2 

(/ 7 'u)^_l_X/ 2 Pr 2 / 2 , 2 + 1/2 (,ddVjfi—i/ 2 Pnh,i—l /2 
“l“Pr 2 / 2 , 2 + 1/2 (0+1 0)^2 “t~ 1^22/2,2—1/2(0 0 — 1 ^ ^2 

-((//«).+! - - {(Hu)i - {hu),w ’’"7~A m) 


(Hu) 


i+1/2 


(Hu)i+i + (Hu)i 


The divergence free condition (|6^ multiplied by l/2pnh,i+i/2 gives 
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Pnh,i+l/2 f f N fzj \ \ Ui -\-Uij^i . . 

' (i/M)i+l - {Eu)i j - - - Pnh,i+l/2 (0+1 - 0) 


{Wi+i + Wi) = 0. (85) 


The sum of relations (|83|), (|84|) and (|85|) gives 

with 


,iPnh- ^(-^^)i+l/2Pnh,i+l/2 (-^^)i—l/2Pnh,j—1/2^ “1“ *^*+1/2 1/2) 


<^j+l/2 — 

di-ll2 = 


Pnh,i^ll2 

2 

Pnh,i—ll2 


'^^i+l/ 2 (^i+l (^i +1 ^ 0 ( 0+1 0 ) j : 


Aa;i_i/2(wi - Wi_i) - {ui - Mi_i)(0 - O-i) 


When the disvergence free condition is satisfied at the boundaries, this proves (iv). 

Assuming the variables are smooth enough, the guantities (ij+i/ 2 ?<^i-i /2 satisfy (ij+ 1 / 2 - 
di-i /2 = 0{Ax^) and we have 


Axi 


^sw,iPnh ((^Hu')i-\--i/2Pnh,i+l/2 (^Huf—-i/2Pnh,i—1/2 j “1“ 0(^Ax') , 


that completes the proof. 


4.3 The fully discrete scheme 

Now we examine the fully discrete scheme that consists in the coupling of the semi¬ 


discrete schemes described in paragraphs 4.1 and 4.2 


4.3.1 Prediction step 


Using the space discretization dehned in paragraph 4.2, we adopt, for the system (53), 
the discretization 

xn+i,2 _ X- - <(Fr+V2- - F"-i/2+), (86) 

where af = Af^/Axi is the ratio between the space and time steps and i^i+i/ 2 ± are given 
by a robust and efficient discretization of the hyperbolic part with the topography. 

For the discretization of the topography source term in the Saint-Venant system, 
several techniques are available. In this paper, we use the hydrostatic reconstruction 
(HR scheme for short) |1], leading to the following expression for the numerical fluxes 


Y'r 

^i+l/2-i^i+l/2+) 


i+1/2- 


rrn 

^i+l/2+ 


rH(xi 

(^i+1/2-’ ^i+l/2+) 

("^i+1/2-’ ^*+1/2+W*+1/2/ 

'^^^■(^i+l/2-> ^i+l/2+) 

-^92: i^i + 1/2-1 ^i+l/2+) 

^‘^^^(^i+1/2-’ ^i+l/2+Wl+l/2/ 


/ 


+ 


9-^ 


y 2 


0 


\ 


(87) 


9Wl 


+ 1 / 2 + 
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where 0 has been used and ^ numerical flux for the Saint-Venant 

system without topography. The reconstructed states 


"^*+1/2- ~ h ^i+i/2+ — ^-^^*+1/2+5i+l/2-|-“i+1 J i 

are defined by 

“ 2;b^i+l/2) + , ^fl!|_i/2+ = {H^+1 + 2:6,i+1 — Zb^i^i/2) + , 

and 

«s,i+i/2 = meix(2^j,2j,i+i). 


( 88 ) 


(89) 

(90) 


4.3.2 Correction step 


For the system (54),(55), we adopt the discretization 


n+l n+1/2 _|_j^ 

Ui = U. ' - J^Vs^,iPnt . 


div^^,i+i /2 (u”+^) = 0 , 


(91) 

(92) 

(93) 


with \/sw,iPnt^ = '^sw,i and div^^,i+i /2 (u^+^) = div^^,i+i /2 

defined by Eqs. (73)-(75). Then, applying dwsw,i+i /2 to (92) and using (93) gives the 

expression for the elliptic equation under the form 


div. 


StD, 2 + 1/2 


^ = ^div,^,i+i/2 (u’^+^/2^ . 


Rn+l 


(94) 


The solution of (94) gives p'^^ and allows to calculate (i7u)+^ using (92). 
Omitting the superscript the expression of 


^sio,fi-|-l/2Pnh diVg^o j^i/2 ( jj ^sw Pnh I : 
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is given by 


^^i+l/2^sw,i+l/2Pnh /\x' Pn/i,i+l/2^ “1“ /\x- Pnh,i—1/2^ 

Pnh,i+3/2 /^ ^ \ Pnh,i+lf2 

“ 

Pnh,i+l/2 / ^ Pnh,i—1/2 j ^ ^ \ 

+ “ ‘■■■' (c- - c-i) 

(Pnh,i+3l2 Pnh,i-\-ll2 Pnh,i+ll2 Pnh,i—l/2\ 

+1 Ayy + Ay ) 

+ (Ch 2 - C,-..) (C,-.. - 0 ) 

Hi+iAxi+i 

Pnh,i+l/2 I Pnh,i+l/2 f^ ^ \2 

+ /yyyyy “ ^‘■> + 

,M(,,._,)(c,-c.-y 

/ AXj+3/2Pnh,i+3/2 + AXj+i/2Pn/i,i+l/2 
* V AXi+iHi+i 

AXi-\.\l2Pnh,i+l/2 ~\~ AXj_i/2Pnh,i—1/2 \ 

Axii/i / ■ 


■ Ax. 


i+l/2 I - 


And it remains to prove the previous relation is consistent with the left hand side of 
Eq. (22). We rewrite Asyj^i^ij 2 Pnh under the form 

Hi+i ( \ Hi / N 

\^nh^i-\-3/2 Pnh,i+ll2J “r \^nh,i-\-l/2 Pnh,i—ll2 


~ ^XiJ^l/2^sw,i-\-l/2Pnh a 

Pnh,i+3/2 f y. r,/- I ^ \ 

1 


^X^ Pnh,i+l/2{Ci+l (i) 


Pnh^i—1/2 ^^ ^ ^ 

h‘+i - 2Ci + C.-1 j 

Pnh,i+3/2 (^ ^ \( /- aA 

zyyAyy h-+' “ ‘■•j 


+ KIlXKj, 


■ + 




(o+i-o)‘ 


AXi+1/2 


/ AXj+3/2Pnh,i+3/2 + AXj+i/2Pn/i,i+l/2 


Axi-i-iETj-i-i 


AXi+i/2Pnft,i+l/2 + AXj_i/2 Pnh,i-l/2\ 

AxiHi ) ’ 

that is indeed a consistent discretization of the left hand side of Eq. (22). 
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In the case of a regular mesh Ax* = Ax = cst, the preceding expression of Asw,i+i/ 2 Pnh 
reduces to 

Ax A^^o Hi^i(^nh,i+3/2 Pn/i,j+l/2^ “1“ ^Pnh,i+l/2 Pnh, 1 — 1 / 2 ^^ 

— Pnh,i+3I2 (^(i+2 “ ‘^Ci+1 + 0 j 

Pnh,i—\l2 ^Ci+1 A 

I Pnh,i+3/2 (^ ^ \ f^ 

\yi+2 — Wl ) Ki+1 ~ ^i) 


H ,+1 

Pnh^i+ll2 


1 


Pnh,i—1/2 I i 

H-77-I Si+l 


H, 


0 0 - 0-1 


Ax' 


2 I Pnh,i+3/2 Pnh,i+l/2 ^ Pnh,i+l/2 Pnh,i—1/2 


H„ 


i+l 


H, 


Remark 4.3 The numerical scheme proposed in this paragraph for the correction step is 
based on a finite difference strategy and hence cannot he extended to unstructured meshes 
in higher dimension. But, based on a variational formulation of the correction step, the 


authors have obtained a finite element version of the scheme (92)-(93) with polynomial 


approximations of the velocities and the pressure satisfying the discrete inf — sup con¬ 
dition. It is not in the scope of this paper to present and evaluate such a discretization 
strategy, it is available in a companion paper m- 

4.3.3 Boundary conditions 

It is difficult to define the boundary conditions for the whole system. Therefore, we first 
impose boundary conditions for the hyperbolic part of the system and then we apply 
suitable boundary conditions for the elliptic equation governing the non-hydrostatic 
pressure Pnh- 


Hyperbolic part The definition and the implementation of the boundary conditions 
used for the hyperbolic part have been presented in various papers of some of the authors. 
The reader can refer to 


Non-hydrostatic part For the non-hydrostatic part, the definition of the boundary 


conditions means to find boundary conditions for Eq. (94) and we adopt the follow 


ing strategy. Notice that other solutions can be investigated since the coupling of the 
boundary conditions between a hyperbolic step and a parabolic step is far from being 
obvious. 

Given flux When the inflow is prescribed - for the hyperbolic part - we impose for the 
elliptic equation (94) a homogeneous Dirichlet type boundary condition. More precisely, 

\n+l/2 


if for i = 0 OT i = I 


1, = Qo is given then we impose = 0- 


= 0. This 
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choice is imposed by the relation (92) in order to ensure on the 

neighbouring cell. 

Given water depth If the water depth is prescribed for the hyperbolic part i.e. for i = 0 


or z = /, ~ ^0 is given then we impose for Eq. (94) a Neumann type boundary 

condition under the form Pnh,1/2 = Pnh,3/2 = 0 or Pnh,i-il2 = Pnh,1+1/2 = 0. 

4.4 The discrete inf — sup condition 

Using a matrix notation for the shallow water gradient operator 

we get 

div,^ = Su"+\ 

where suitable boundary conditions are assumed. Therefore, defining 

A = diag(i/"^), 

the fully discrete scheme obtained from (|86j),(91 )-(93) can be rewritten under the form 


'it 0 0 

0 ii 

0 BA 0 


pjn+l \ / H 


[Hu 


.n+l 






(95) 


where Dh,Dhu refer to the numerical discretization of the hyperbolic part. The fact 
that the matrix BAB^ is invertible is related to the inf-sup condition [15], see [1] where 
this property is investigated in details. 


Remark 4.4 Instead of the scheme (95), a fully implicit version - including the hyper¬ 


bolic part - may he considered. But such a discretization would imply to have an implicit 
treatment of the hyperbolic part of the proposed model corresponding to the Saint- Venant 
system. And an efficient and robust implicit solver for the Saint-Venant system is hardly 
accessible. 


4.5 Stability of the scheme 


For the numerical scheme detailed in paragraphs |4.3.1 and 4.3.2 we have the following 
proprosition. 


Proposition 4.5 The scheme (86),(91 )-(93) 


(i) preserves the nonnegativity of the water depth > 0, Vz, Vzz, 

(ii) preserves the steady state of the lake at rest, 

(in) is consistent with the model (§-(§. 
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Proof of prop. 4.5| (i) The statement that preserves the nonnegativity of the water 
depth means exaetly that 


^H 0, Hi^i , Mj-i-i) ^u(^H.i—1 , 'Uj—i, Hi 0, uf) ^ 0, 

for all ehoices of the other arguments. From ([8^, (87), we need to check that 

- 1 / 2 +) 




whenever Hf = 0. And this property holds since from (89),(90) Hi = 0 implies Hi+i/ 2 - = 


H 


i-l/2+ 


= 0 . 


(a) When u2 = 0 for all i, the properties of the hydrostatic reconstruction technique 
ensures 

rpn _ pn 

■^ 2 + 1 / 2 - “ -^2-1/2+5 


in (86) and hence = Xf moreover the scheme (91),(92),(94) gives 

^n+l _ ji^n+ 1/2 

proving the scheme is well-balanced. 

(Hi) The numerical flux X being consistent with the homogeneous Saint-Venant sys¬ 
tem, the hydrostatic reconstruction associated with X gives a consistent discretization of 


the Saint-Venant system with the topography source term. The discretizations (92),(93) 
being obviously consistent with the remaining part, this proves the result. ■ 

4.6 

When H tends to 0, the correction step (92) is no longer valid and we propose a modified 
version of (92),(93) under the form 

1 


+1 = 


diysui,i+l/2 


u 


'”+^) = 0 , 


and 


Xxi +1 


vye I = Ti ’^+1 — 71”+^ 

^sw,it'nh li rnh,i+ll2 t'nh,i-l/2 


^71+1 ''sw 


V" ■p^t^\ = 

^sw^it'nh I 2 


' TTn-\-l I -^77/1,2+1/2 \^^2+l J ' T^nh,i—1/2 >2—1 J 

■^ 2,6 V 

~ prn+l {)^^i+^/^Pnh,i+l/2 + ^^i-^/'^Pnh,i-l/2^ ’ 

P 


with e: being a constant £ = cst > 0 and H^ = max(i7, e:). 
In order to ensure the total pressure 

Pu I - 

-H + 
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remains non negative, we add the constraint 

9 


pSI+ 1/2 = 0 I -ff.?/) + P"j+i/2 < ». 


to the solution of the elliptic equation (94). Notice that in all the numerical tests 


presented in this paper this constraint is not active meaning the total pressure remains 
non negative. 


5 Fully discrete entropy inequality 


We have precised in paragraph |4.3| a general scheme for the resolution of the non¬ 
hydrostatic model. In this paragraph we study the properties of the proposed scheme in 
the context of one particular solver for the hyperbolic part, namely the kinetic solver, 
since it allows to ensure stability properties among which are entropy inequalities (semi¬ 
discrete and fully discrete) |5]. 

Looking for a kinetic interpretation of the HR scheme, we would like to write down 


a kinetic scheme for Eq. (45) such that the associated macroscopic scheme is exactly 


(86)-(87) with the dehnitions (88)-(90). 


We drop the superscript "■ and keep superscripts n-|- 1 and u-|- 1/2. We denote M* = 
M{Hi,Ui,^), Mj+i/2- = ^*+1/2+ = dL(iLj+i/2+, Mj+l, 0; fi~^ ^ — 

(0) we consider the scheme 


^ — Mi — ai ^^lg<oMj_|_i/2+ + ^lg>o^i-i-i/2- + ^Mi+i/2- 

~^^^>oMi-i/2- — ^l^<od4i_i/2+ — (5Mj_i/2+^ • 


(96) 


In this formula, (5Mj_|_i/2± are defined by 

SMi^i/2- = {^ — Ui){Mi — Mj+i/2-), 5 M(+i/2+ = {^ — Wi+i)(Mj+i — Mj_|_i/2+), 

and are assumed to satisfy the moment relations 

t2 


/ SMi+ii2- = 0 , 

Jr 

f dMi_i/2+ = 0 , 

Jr 

Dehning the update as 


f SMi+y2- = 9^ - 9 1 

? de = g:S _ 


(97) 

(98) 


\ 71+1/2 ^ 


Hu 


Ax. 


*i+l/2 


and 




Ax 


^ Q) ,x,i)dxdi, 

(99) 

’*+1/2 


/(r+ ! -,x,^)dx, 

(100) 
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the formula (99) can then be written 


H 

Hu 


n+l/2 




and using (71) we also dehne 

= H,w, - a,(n;,+i/2F^;r+i/2 “ 


( 101 ) 


( 102 ) 


with 


Fh^+I/2 — / ^(ls>oMi+l/2- + l^<0^i+l/2+)<^^- 


Finally, relations (96)-(102) give an explicit formula for the prediction step (86) and we 
have the following proposition that is proved in |ni Corollary 3.8] (the two constants Vm 
and Cjs are precised therein). 

Proposition 5.1 Let Vm and Cp he two constants. For ai small enough, the numerical 
scheme (86) based on the kinetic description (96)-(102) and the HR technique (89),(90) 
satisfies the fully discrete entropy inequality 

(103) 


“1“ C^(cyiVyn^ I l) ) 5 


with Gi+l /2 = Ghyd,i+l/2 + FH,i+l/2- 


i+1/2 


and 


G 


hyd^i-\-ll2 


'C<0 


fH(Mi^i/2+, Zh^i^i/ 2 ) d^ + / ^H(Mi^i/2-, Zb^i^i/ 2 ) d^. (104) 


'^>0 


Remark 5.2 Let us notice that the quadratic error term has the following key properties: 
it vanishes identically when z = cst (no topography) or when cxj —?• 0 (semi-discrete 
limit), and as soon as the topography is Lipschitz continuous, it tends to zero strongly 
when the grid size tends to 0 (consistency with the continuous entropy inequality Qt 
even for non smooth solutions. 


For the correction step, we have the following discrete energy balance. 


Proposition 5.3 The numerical scheme (92),(93) satisfies the following inequality 

v(xr‘) < f,(xp''^) 



Cj((Ai”)" + C2(Aiy). 


Moreover, for cxj ~ 
have 


1, At"' small enough and assuming suitable boundary conditions, we 


Yi Aii (i)(A'"+') - MxYL) < 0. 


26 















Corollary 5.4 The numerical scheme detailed in paragraphs 4-3.1 and \4.3.'^ satisfies 
the fully discrete entropy inequality 

< fj{Xi) — CTi^Gj+i/2 — 

+Cy{aiVmY (^g{zb^i+i-Zb,iY+g{zb,i-Zb^i_iY'^ (af+C2Ar), 
with Gi^i /2 = GiJ^i /2 + {Hu)i^i/ 2 Pnh,i+l/ 2 - 


Proof of prop. 5.1 With rjhydiX) defined by (10); we start from the inequality 
Vhydi^^i ) ^ Vhydi^^i) G}iyd^i—\l2 


“1“ GyiyJiV'ffi) “1“ ^ : 


(105) 


that is proved in Corollary 3.8j. Equation (105) corresponds to a fully discrete entropy 
inequality for the Saint- Venant system including the topography source term. 


Multiplying (102) by Wi leads to 




with 


Tjn+l/2 jj 


w- 


+ - H.) - 2 


h: 


n+l/2 


/ n+l/2 \2 

K -Wi) , 


and 


wfi wfi 

n,, ipkin TTifcm _ i+1/2 j^kin i—1/2 j^kin 

Wi+l/2Wir fj^i^lj2 Wi_i/2Wir fji_i^2 — 9 ^ H,i+l/2 9 ^ H,i-l/2 


^*+1/2 \ rpkin Wi-i/2 


+ W.+1/2F^:?+1/2 ) - W,.y2FrU/2 - 2 

The sum of the two previous relations gives 


h: 


n+l/2 


2 xn+l /2 Hi 2 , ^ f'^i+1/2 Tpkin rpkin 


I 2^’^+i/^ -‘-‘i 2 I 


JpUlTL ^ T?' 

^ ^ H,i+\I2 ^ -^H,i-l/2 


72x74+1/2 

< T K?+1/2]_ K+l - - ? K?-./J+ (K-i - K'-l)" + - Vlif. 


( 106 ) 


It remains to estimate, when ^ quantity 


H 


n+l/2 


/ n+l/2 n2 

[Wi -Wi) , 
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in the r.h.s. of Eq. (106). 
We have 

n+1/2 

Wi -Wi = 


n+1/2 


h: 


- iHw)i - - H,)w, 


h: 


.„+l/2 '^i+ 1 / 2 ) Wi_i/2)\, 


and hence 

Tjn+1/2 

^L_f„“+V2 




h: 


Therefore, for ai small enough, the r.h.s. of Eq. (106) is non positive with 

/7/?^ 7/?^ 

'^*..,,2 I „ ( *+1/2 jpkin i—ll2 jpkin '\ ^ /^2„ /-\ r^2 „ \ 

^- Wi) - Y^i -^^’*+1/2- lf—PH,i-l/ 2 ) < -CiCTi(l -6211*). 


T/ie previous relation coupled with (105) gives the result. 


Proof of prop. 5.3| We start from relation (92) multiplied by , this gives 

ri+1/2 ^n+1/2 


H, 




n+1 


H, 


*u? 


2 

/ 5 


with the notation = u.u. 

Omitting in this part the superscript and as in the the proof of prop. 4.2, simple 
manipulations give 


with 

di+ 1/2 = 
di-\i2 = 


,iPnh- (^(^Hu)i-^-i/2Pnh,i+l/2 l/2Pn/i,i—1/2^ “1“ c(i+l/2 1/2; 

u,+,-Ui ^ + 2zk,i))), 


Pnh.,i-\-ll2 

2 

Pnh,i—ll2 


AXi+i/2{Wi+i - Wi) - 
AXi_i/2{Wi -Wi_i) - 


^i—1 


{Hi + 2zb^i — {Hi_i + 2zb,i-i))\ 


2 t xjt t xj 2 

proving the result. 

Assuming the variables are smooth enough, the quantities (ij+i/ 2 ;<^i-i /2 satisfy (ij+ 1/2 
di-i /2 = 0{Ax^) and we have 


Axi 


■^sw,iPnh (^(.Hu)i^i^2Pnh,i+l/2 (.Hu)i—i/2Pnh,i—l/2'j “1“ C^(Ax) , 


that completes the proof. 

Notice that in the limit AW —)■ 0, cTj ~ 1, the inequality 


Tjn+l/2 2 

*U»+1 _ u”+V2\ 


J + o'i(di^i/2 — di- 1 / 2 ) < 0, 

holds, meaning the correction step ensures a decrease of the entropy. ■ 

Proof of corollary |5.4| The sum of the two inequalities obtained in props. 5A and T_3 
gives the result. ■ 
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6 Analytical solutions 


Stationary and time dependent analytical solutions are available for the model (18)-(20), 


see m and references therein. In this section we only briefly recall some of them, they 
will be very useful to evaluate the properties of the proposed numerical scheme, see 
paragraph 


6.1 Time dependent analytical solution 


6.1.1 Parabolic bowl 



The functions defined by 



H{x,t) = 

max ^Hq “ Y f{ti)dti^ , 0 j , 

(107) 

u{x,t) = 

f{t)^H>0, 

(108) 

w{x,t) = 

b2Xf{t)lH>0, 

(109) 

Zb{x) = 

bi + 

( 110 ) 

Pnhix,t) = 

2 

( 111 ) 

s{x,z,t) = 


(112) 


where Hq > 0 , 61,62 are constants and the function / satisfies the ODE 

fO TO 


^ + b2{g + 62/^) f{ti)dti = 0, /(to) = f, e 


are solutions of the system 

A 

dx 




dt 

d 


d 


^(Hu) + ^{Hu^ + ) = -{gH + 


dt dx 

— (/Dh) + —{Hum) = 2p^f^ + Hs, 


dzh 


dx 


dt 


dx 


am) + 


dx 


dx 


(113) 

(114) 

(115) 

(116) 
(117) 


tion (116) 


that corresponds to (18)-(19) completed with a source term Hs in the momentum equa¬ 


6.1.2 Solitary wave solutions 


The system (|18|)-(|19|) admits solitary waves having the form 

' X — Cot 


H = Ho + a 


(.ech(i 


I 


(118) 
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“ = Co (^1 - ^ , 

_ acod , /X — cot\ ^ , f X — Cot 

w = ——— seen I --- seen I --- 

l±l I J I 


+i^seeh 


X Cot ^ ^ ff (^ ^ot 


I 


seeh 


where ip' denotes the derivative of funetion (p, 


(119) 

( 120 ) 


( 121 ) 


Co = 


I gHi 


d \l P-Hi' 


a = 


-^0 


P-Hl' 


and {d, I, Ho) G are given eonstants with I > Ho > 0. 


7 Numerical simulations 


A eomplete validation of the proposed nunierieal teehnique is not in the seope of this pa¬ 
per and will be investigated in a fortheoniing paper. We foeus on two typieal situations, 
deseribed in paragraphs 6.1.1 and 6.1.2| where analytieal solutions exist. 

For the nunierieal test, we use a 2”'^ order extension of the spaee diseretization for 
the predietion step. The seeond order extension is built as in |7]. For the seeond order 
extension of the time seheme, we use a Heun type seheme, see m- 


7.1 The parabolic bowl 


(see also m) is 


At the discrete level, the analytical solution given in paragraph 6.1.1 
particularly difficult to capture. Indeed, it is a non stationary solution and the flow 
exhibits wet/dry interfaces all along the simulation. 

With the parameter values Ho = 1, a = 1, bi = 0, b 2 = 1, P = 0 over the geometrical 
domain [—2,2] and with the initial conditions (see Fig. 


r m dt = 

Jio V gb2 

fM = /° = 0, 


we have calculated - with a simple Runge-Kutta scheme - the solution of the ODE (113). 


The solution has been calculated with a very hne time discretization and thus can be 


considered as a reference solution, very close to the analytical solution of Eq. (113). This 


means we have at our disposal an analytical solution for the system (114)-(117). 
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Figure 2: Initial conditions for the simulation of the “parabolic bowf’ (parabolic bottom, 
water depth and null horizontal velocity). 


To illustrate the behavior of the solution in such a situation, we give over Fig the 
variations along time of the water depth at xq = 0.8 m for a mesh of 80 cells which is a 
rather coarse mesh. 

In order to evaluate the convergence rate of the simulated solution towards the an¬ 
alytical one, we plot the error rate versus the space discretization. We have plotted in 
Fig. I^the log(L^ — error) over the water depth at time T = 10 seconds - corresponding 
to more than 5 periods - versus logQiQ/hi) for the first and second-order schemes and 
they are compared to the theoretical order. We denote by hi the average cell length, 
ho is the average cell length of the coarser space discretization. These errors have been 
computed on 6 meshes with 20, 40, 80, 120, 300 and 400 cells. 

For this test case, the errors due to the time and space schemes are combined so the 
convergence rate of the simulated solution towards the analytical one is more difficult to 
analyze. Despite this fact, it appears that the computed convergence rates are close to 
the theoretical ones. It emphasizes the performance of the proposed numerical technique. 

Over Fig|^ we have plotted the variations along time of the log(L^ — error) over the 
water depth for the mesh with 80 cells calculated at node Xq = 0.8 m i.e. the quantity 

Hsimip^Oit) (Xq, f) I \ 

Hsim{xQii) ) 

It appears that this quantity does not increase with time and remains bounded. 


f log ( 100 


7.2 The solitary wave 


As in the previous paragraph, we are able to examine the convergence rate of the simu¬ 
lated solution towards the analytical one (given in paragraph 6.1.2). 
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Figure 3: Parabolic bowl: variations of t h-)■ if (0.8, t) - analytical solution and simulated 
one with the first order (space and time) and second order (space and time) schemes. 



Figure 4; Parabolic bowl; convergence rates to the reference solution, order schemes 
(space and time) and 2”'^ order schemes (space and time), theoretical order. 
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Time (s) 


Figure 5; Parabolic bowl; variations along time of the error at node xo = 0.8 m, P* 
order schemes (space and time) and 2"''^ order schemes (space and time). 


We consider the analytical solution corresponding to the choices Ho = 1 m, / = 1.7 
m, d = 1 m. These choices lead to a = 0.5291 m and Cq = 3.873 m.s“^. We compare the 
analytical solution and its simulated version at time t = 6 s. 

We have plotted in Fig. [^the log(L^ — error) over the water depth at time T = 6 
seconds versus log{ho/hi) for the first and second-order scheme and they are compared 
to the theoretical order. These errors have been computed on 7 meshes with 80, 120, 
200, 400, 800, 1600 and 3200 cells. 


For the curves obtained over Fig. the soliton is - at the initial instant - in the 
fluid domain meaning the numerical treatment of the boundary condition does not play 
a crucial role. Figure is similar to Fig. except that the soliton is not within the 
fluid domain at the initial instant but enters the channel by the left boundary. The 
convergence order of the scheme is examined at time T = 10 seconds and the soliton 
arrives within the fluid domain after 4 seconds of simulation. Hence the soliton propa¬ 
gates during 6 seconds within the domain corresponding to the same situation as Fig. 
Following 4.3.3 we have imposed a given flux at the entry of the domain. We notice 
that the convergence orders obtained over Figs. and are similar. 


8 Conclusion 

In this paper we have proposed a robust and efficient numerical scheme for a non¬ 
hydrostatic shallow water type model approximating the incompressible Euler system 
with free surface. 
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Figure 6; Soliton (interior of the domain); convergence rates to the reference solution, 1®* 
order schemes (space and time) and 2”*^ order schemes (space and time), theoretical 
order. 



0.5 r 


0- 


-0.5- 

E 


o 

c 

-1 

_l’~ 


c 

-1.5- 

o 


0 

o 

-2 

O) 

o 

-2.5- 


-3 



2 2.5 3 3.5 

log (h^/h.) 


4 


Figure 7: Soliton (entering in the domain); convergence rates to the reference solu¬ 
tion, P* order schemes (space and time) and 2"''^ order schemes (space and time), 
theoretical order. 
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The correction step is the key point of the scheme and especially the discretization 
of the shallow water type version of the gradient and divergence operators. A finite 
difference strategy has been used to discretize this correction step. In order to be able 
to treat 2d flows on unstructured meshes, a variational approximation of the correction 
step is required, it is presented and its numerical performance is evaluated in [Tj. 
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